关键词：两相流;VOF方法;界面重构
Advances in research of VOF method
ZHANG Jian1, FANG Jie1, FAN Boqin2
(1.College of Water Conservancy and Hydropower Engineering, Hohai Univ.,
Nanjing 210098, China;
2. Water Resources and Hydropower Technology Consulting Center of
Zhejiang Province, Hangzhou 310012, China)
Nanjing 210098, China;
2. Water Resources and Hydropower Technology Consulting Center of
Zhejiang Province, Hangzhou 310012, China)
Abstract: An introduction was given to the principle of the VOF meth od in the simulation of interface variation in twophase flows.Some works for improvement of the calculation precision and the effect of interface reconst ruction of the VOF method and some research results were summarized. Finally, it was pointed out that the VOF method could be widely applied to numerical simula tion of multidimensional flow fields if some theoretical improvements were made.
Key words: twophase flow; VOF method;interface reconstruction
Key words: twophase flow; VOF method;interface reconstruction
1VOF方法及原理
VOF方法的基本原理是通过研究网格单元中流体和网格体积比函数F来确定自由面,追踪流体的变化,而非追踪自由液面上质点的运动.VOF方法可以处理自由面重入等强非线性现象,所需计算时间短、存储量少,但在处理F的变化时稍显繁琐,有一定人为因素.VOF方法根据体积比函数F来构造和追踪自由面.若F=1,则说明该单元全部为指定相流体所占据；若F=0,则该单元为无指定相流体单元；当0＜F＜1时,则该单元称为交界面单元.假定流场中任意点(x, y),定义函数f(x,y,t)如下：
采用图1所示的交错网格,阴影表示液体部分,时间采用一阶差分格式,方程(2)差分形式为
量(图1箭头方向表示流通量方向),其积分形式为
2 VOF方法研究进展
自Hirt和Nichols提出VOF方法以来,为了模拟强非线性流动,国内外学者从方程差分格式和自由面流体传输两方面着手,结合自由面重构格式对VOF方法进行了联合改进.Noh等［2］提出的SLICVOF方法采用分裂算法,规定界面法向与坐标轴平行,人为性较强.Hirt和Nichols采用施主受主格式:图1中(i,j)单元右边界的流通量由边界相邻单元的体积比函数F以及边界速度共同确定,具体见式(5),自由面方向由相邻8个单元的F值规定.Cook等［3］在此基础上推广至三维,对输运方程在空间上离散采用二阶迎风格式,网格划分采用非正交自适应网格,流体输运规则仍采用施主-受主法则.
Puckett等［4］为保证质量守恒,对输运方程采用方向分裂算法,将方程(2)改写成方程(6)和(7)，再沿两个方向交替进行数值计算.该方法对自由面则用离散的线段来表示,通过相邻网格的体积分数来预测自由面的法向,近似重构得到的结果相对较粗糙.Gao等［5］通过计算认为在一个时步内,交替方向重复循环计算可得到更的F值.
FCT-VOF方法由Boris和Book［6］提出并经Zalesak［7］改进扩展至多维.扩展方法有两种：①基于数值流通量守恒来求解体积函数,不分裂算子.二维x方向上右边界的流体体积量由式(8)和式(9)求解：
其中低阶流通量采用顺风格式，高阶流通量采用逆风格式,
为流量校正因子.该格式既可消除逆风格式的耗散效应,又可抑制顺风格式的不稳定.②采用方向分裂算法.Rudman［8］通过验证发现后者模拟效果更好，但缺点是在进行交替方向计算时得到的中间值F可能会大于1.0.由于在进行另一方向的计算时要用到F,因此原来满流体单元内部的值会小于1.0,产生一些非物理的破碎区域.Rudman［8，9］在此基础上针对交替方向计算可能会出现F大于1.0的情况,对每个方向计算的是流体体积,再利用Lagrange的观点修正网格体积,两者比值即为中间体积比函数F,对每个方向进行计算,zui后得到网格下一时刻的流体体积比函数F.此外Rudman还对多相流体交界面建立数学模型,认为对于多相流体交界面均可采用流体输运方程,只是其中流体体积比函数F的定义不同.
Youngs［10］提出的PLIC方法对于满单元采用一阶顺风格式.对于自由面单元,先计算自由面法向,*次引入自由面直线与x轴夹角β,根据β判断自由面类型,并由此计算流过四周边界的流体体积,来修改网格内的流体体积,考虑了相邻界面流体的细微输运,计算结果非常精细.重构方法则注意了自由面直线与四边的交点.Guignard［11］通过在自由面的斜线段上设置标记点,根据标记点的移动方位来判断得到下时刻的F,与PLIC方法相比,这样做的优点是可以采用更大的时间步长.1991年Liang［12］引入对角法向单元,从而将VOF法从正交坐标系拓展到非正交坐标系.以上算法中自由面单元的连接基本上都是协调的,多数采用一阶精度的施主受主模型.
Ashgriz和Poo［13］提出的FLAIR-VOF方法也采用分裂算法.对于相邻单元为满单元,输运方程采用一阶迎风格式;对于含自由面单元,则根据相邻网格的体积分数求出自由面斜线段,再根据斜线段斜率以及网格边界速度方向来决定边界上的流通量,输运方程在空间离散可以达到二阶精度,修正后的Lax-Wendoff格式仅适用于均匀网格.王志东［14］对该方法进行了修改,使其适用于非均匀网格.李谊乐等［15］采用TVD差分格式,根据Hirt-VOF方法原则采用1／4和3／4单元表示自由面网格中流体的形状,使自由面的TVD格式可以达到二阶精度，即可以消除自由面剧烈变化时Lax-Wendr off格式引起的震荡和一阶迎风格式的耗散效应.Lafaurie等［16］将三维输运方程离散,根据自由面的法向和边界面流速方向夹角值来确定边界面的流通量.夹角值由边界面相邻单元的流体体积确定，显然这样规定较简单.Gueyffier等［17］将PLIC方法推广至三维,边界面上流通量由相邻网格单元流体体积数及边界面上中心速度和相邻单元中将要流过边界的质点速度共同确定;界面构造则是寻找运动界面所穿过的网格范围,仅进行与运动界面相交及相邻网格的运算,界面的位置和方向采用几何和代数方法迭代得到,并推导出了自由面方程的近似解析解.
Ubbink和Issa［18］根据Gaskell和Lau 的可压缩性构造思想,提出了CICSAM方法.控制体内流体体积为连续流构造的流通量和可压缩性思想构造的流通量的平均值,利用精细网格内界面的迁移方法得到界面在网格内的运动位置和法向.该方法采用半隐格式处理,根据式(10)和(11),不用分裂算子就可扩展至多维.方程采用C-N格式离散,并对网格边界上的流通量进行修正,既保证有界,又维持体积函数守恒,并且适用于无结构网格.
式中：P，K，I，VP，FI分别代表控制体中心位置、边界数目、控制体边界面位置、控制体体积、第I边界上的流通量.βI通过NVD方法确定.
Dendy等［19］对CICSAM方法做了两点改进：*，输运方程采用显格式,修正流体输运公式,采用多维连续流格式代替CICSAM方法中的ULTIMATE-QUICKEST高阶一维格式,适用于四边形、三角形网格.对于连续水流区域的控制体K,第I面流通量δFI表示如式(12).第二，继续使用CICSAM方法的C-N项,但采用Thuburn流量因子,对于包含其他相流体的不连续流体区域,则用高阶迎风显格式和光滑连续隐格式相结合代替C-N格式,以便在大库朗特数情况下也可以计算,且适用于无结构网格.
式中：VK，AI分别为第K控制体体积和第I边界面面积.
Hong［20］提出了一种斜线段界面重构法,采用水平相邻3个单元的F值来确定斜率,方程系数由相邻8个单元的F值来确定.Kim等［21］采用二次曲线来构造自由面,做得更为精细.图2为几种常用VOF法对实际流体进行构造,其中图2(a)为实际交界面.分别采用Flair,Youngs,Hirt方法的构造规则对图2（a）进行构造，结果如图2（b）、图2(c)和图2（d）所示.图2(b)和图2（c）采用斜线段来构造自由面,得到的结果与实际情况更吻合.
自Hirt和Nichols提出VOF方法以来,为了模拟强非线性流动,国内外学者从方程差分格式和自由面流体传输两方面着手,结合自由面重构格式对VOF方法进行了联合改进.Noh等［2］提出的SLICVOF方法采用分裂算法,规定界面法向与坐标轴平行,人为性较强.Hirt和Nichols采用施主受主格式:图1中(i,j)单元右边界的流通量由边界相邻单元的体积比函数F以及边界速度共同确定,具体见式(5),自由面方向由相邻8个单元的F值规定.Cook等［3］在此基础上推广至三维,对输运方程在空间上离散采用二阶迎风格式,网格划分采用非正交自适应网格,流体输运规则仍采用施主-受主法则.
Puckett等［4］为保证质量守恒,对输运方程采用方向分裂算法,将方程(2)改写成方程(6)和(7)，再沿两个方向交替进行数值计算.该方法对自由面则用离散的线段来表示,通过相邻网格的体积分数来预测自由面的法向,近似重构得到的结果相对较粗糙.Gao等［5］通过计算认为在一个时步内,交替方向重复循环计算可得到更的F值.
FCT-VOF方法由Boris和Book［6］提出并经Zalesak［7］改进扩展至多维.扩展方法有两种：①基于数值流通量守恒来求解体积函数,不分裂算子.二维x方向上右边界的流体体积量由式(8)和式(9)求解：
其中低阶流通量采用顺风格式，高阶流通量采用逆风格式,
为流量校正因子.该格式既可消除逆风格式的耗散效应,又可抑制顺风格式的不稳定.②采用方向分裂算法.Rudman［8］通过验证发现后者模拟效果更好，但缺点是在进行交替方向计算时得到的中间值F可能会大于1.0.由于在进行另一方向的计算时要用到F,因此原来满流体单元内部的值会小于1.0,产生一些非物理的破碎区域.Rudman［8，9］在此基础上针对交替方向计算可能会出现F大于1.0的情况,对每个方向计算的是流体体积,再利用Lagrange的观点修正网格体积,两者比值即为中间体积比函数F,对每个方向进行计算,zui后得到网格下一时刻的流体体积比函数F.此外Rudman还对多相流体交界面建立数学模型,认为对于多相流体交界面均可采用流体输运方程,只是其中流体体积比函数F的定义不同.Youngs［10］提出的PLIC方法对于满单元采用一阶顺风格式.对于自由面单元,先计算自由面法向,*次引入自由面直线与x轴夹角β,根据β判断自由面类型,并由此计算流过四周边界的流体体积,来修改网格内的流体体积,考虑了相邻界面流体的细微输运,计算结果非常精细.重构方法则注意了自由面直线与四边的交点.Guignard［11］通过在自由面的斜线段上设置标记点,根据标记点的移动方位来判断得到下时刻的F,与PLIC方法相比,这样做的优点是可以采用更大的时间步长.1991年Liang［12］引入对角法向单元,从而将VOF法从正交坐标系拓展到非正交坐标系.以上算法中自由面单元的连接基本上都是协调的,多数采用一阶精度的施主受主模型.
Ashgriz和Poo［13］提出的FLAIR-VOF方法也采用分裂算法.对于相邻单元为满单元,输运方程采用一阶迎风格式;对于含自由面单元,则根据相邻网格的体积分数求出自由面斜线段,再根据斜线段斜率以及网格边界速度方向来决定边界上的流通量,输运方程在空间离散可以达到二阶精度,修正后的Lax-Wendoff格式仅适用于均匀网格.王志东［14］对该方法进行了修改,使其适用于非均匀网格.李谊乐等［15］采用TVD差分格式,根据Hirt-VOF方法原则采用1／4和3／4单元表示自由面网格中流体的形状,使自由面的TVD格式可以达到二阶精度，即可以消除自由面剧烈变化时Lax-Wendr off格式引起的震荡和一阶迎风格式的耗散效应.Lafaurie等［16］将三维输运方程离散,根据自由面的法向和边界面流速方向夹角值来确定边界面的流通量.夹角值由边界面相邻单元的流体体积确定，显然这样规定较简单.Gueyffier等［17］将PLIC方法推广至三维,边界面上流通量由相邻网格单元流体体积数及边界面上中心速度和相邻单元中将要流过边界的质点速度共同确定;界面构造则是寻找运动界面所穿过的网格范围,仅进行与运动界面相交及相邻网格的运算,界面的位置和方向采用几何和代数方法迭代得到,并推导出了自由面方程的近似解析解.
Ubbink和Issa［18］根据Gaskell和Lau 的可压缩性构造思想,提出了CICSAM方法.控制体内流体体积为连续流构造的流通量和可压缩性思想构造的流通量的平均值,利用精细网格内界面的迁移方法得到界面在网格内的运动位置和法向.该方法采用半隐格式处理,根据式(10)和(11),不用分裂算子就可扩展至多维.方程采用C-N格式离散,并对网格边界上的流通量进行修正,既保证有界,又维持体积函数守恒,并且适用于无结构网格.
式中：P，K，I，VP，FI分别代表控制体中心位置、边界数目、控制体边界面位置、控制体体积、第I边界上的流通量.βI通过NVD方法确定.
Dendy等［19］对CICSAM方法做了两点改进：*，输运方程采用显格式,修正流体输运公式,采用多维连续流格式代替CICSAM方法中的ULTIMATE-QUICKEST高阶一维格式,适用于四边形、三角形网格.对于连续水流区域的控制体K,第I面流通量δFI表示如式(12).第二，继续使用CICSAM方法的C-N项,但采用Thuburn流量因子,对于包含其他相流体的不连续流体区域,则用高阶迎风显格式和光滑连续隐格式相结合代替C-N格式,以便在大库朗特数情况下也可以计算,且适用于无结构网格.
式中：VK，AI分别为第K控制体体积和第I边界面面积.
Hong［20］提出了一种斜线段界面重构法,采用水平相邻3个单元的F值来确定斜率,方程系数由相邻8个单元的F值来确定.Kim等［21］采用二次曲线来构造自由面,做得更为精细.图2为几种常用VOF法对实际流体进行构造,其中图2(a)为实际交界面.分别采用Flair,Youngs,Hirt方法的构造规则对图2（a）进行构造，结果如图2（b）、图2(c)和图2（d）所示.图2(b)和图2（c）采用斜线段来构造自由面,得到的结果与实际情况更吻合.
3VOF法的应用和前景
VOF方法简单而有效,在实际工程应用中取得了满意的计算结果.李然等［22］利用三维VOF模型模拟水利工程中天然河流水面线,计算分析得到的洪水对泄洪渠的影响程度和范围与实测结果相比误差很小;姚朝辉等［23］利用HIRT-VOF方法对核电站热力管系统中存在的蒸汽和冷水混流中出现的蒸汽泡冷凝溃灭过程进行数值模拟,分析了热力管系统中出现蒸汽泡和冷水交汇时气泡溃灭的变化过程,指出热力管中蒸汽泡溃灭将产生巨大的水锤压力；包光伟［24］采用HIRTVOF方法对火箭燃料箱内的液体在关机时段的大幅晃动响应过程进行了模拟,结果表明该方法计算液体飞溅和破碎现象的能力较强;李梅娥等［25］通过SOLA-VOF方法对金属铸造成型过程中的流体运动规律进行数值模拟,得出改变浇铸系统浇口和改变液体温度应力分布可以提高液面稳定、改善充型形态的结论;戴会超等［26］采用HIRT-VOF方法较准确地模拟出淹没水跃的大尺度紊流结构及其微观结构;Hur DongSoo等［27］采用三维SOLA-VOF方法对无碎波海域中不对称建筑物周围的水流流态进行数值模拟,准确俘获到建筑物附近的海域波型;Hieu Phung Dang等［28］采用HF-VOF方法对斜坡海岸上的波浪翻滚现象进行模拟,指出采用N-S双方程耦合求解可改进自由面单元中的压力和速度插值;Tseng等［29］采用PLIC-VOF方法模拟微观管网系统中微观井的充液过程,指出井的形状和液体的表面张力对充液结果影响很大.
近年来,随着我国对西部水电的开发,输水管线越来越长,布置形式日趋复杂,由于长时间停机造成管道中聚集大量气体,短时间无法顺利排出,在机组启动时,可能产生比常规水锤大得多的压力;此外,在一些大型水电站输水系统的方案比选中,提出了一种新的布置形式——变顶高尾水洞方案,当输水系统发生水力瞬变时,变顶高尾水洞尾水出流为明满交替流;采用无结构网格进行数值模拟分析,得到的结果将会更符合实际;另外,随着城市建设不断发展,供排水系统日趋复杂,当突降暴雨时,地面积水由于排水系统中存在气体而无法顺畅排出,水流流态由于地下系统气、液交界面的存在变得相当复杂,甚至会对地下排水系统的安全构成威胁.由此可见,很多涉及水、气两相交界,互相推动的水力瞬变问题亟待解决,而VOF方法作为处理自由面边界条件的重要方法,将在未来发挥越来越大的作用.
VOF方法简单而有效,在实际工程应用中取得了满意的计算结果.李然等［22］利用三维VOF模型模拟水利工程中天然河流水面线,计算分析得到的洪水对泄洪渠的影响程度和范围与实测结果相比误差很小;姚朝辉等［23］利用HIRT-VOF方法对核电站热力管系统中存在的蒸汽和冷水混流中出现的蒸汽泡冷凝溃灭过程进行数值模拟,分析了热力管系统中出现蒸汽泡和冷水交汇时气泡溃灭的变化过程,指出热力管中蒸汽泡溃灭将产生巨大的水锤压力；包光伟［24］采用HIRTVOF方法对火箭燃料箱内的液体在关机时段的大幅晃动响应过程进行了模拟,结果表明该方法计算液体飞溅和破碎现象的能力较强;李梅娥等［25］通过SOLA-VOF方法对金属铸造成型过程中的流体运动规律进行数值模拟,得出改变浇铸系统浇口和改变液体温度应力分布可以提高液面稳定、改善充型形态的结论;戴会超等［26］采用HIRT-VOF方法较准确地模拟出淹没水跃的大尺度紊流结构及其微观结构;Hur DongSoo等［27］采用三维SOLA-VOF方法对无碎波海域中不对称建筑物周围的水流流态进行数值模拟,准确俘获到建筑物附近的海域波型;Hieu Phung Dang等［28］采用HF-VOF方法对斜坡海岸上的波浪翻滚现象进行模拟,指出采用N-S双方程耦合求解可改进自由面单元中的压力和速度插值;Tseng等［29］采用PLIC-VOF方法模拟微观管网系统中微观井的充液过程,指出井的形状和液体的表面张力对充液结果影响很大.
近年来,随着我国对西部水电的开发,输水管线越来越长,布置形式日趋复杂,由于长时间停机造成管道中聚集大量气体,短时间无法顺利排出,在机组启动时,可能产生比常规水锤大得多的压力;此外,在一些大型水电站输水系统的方案比选中,提出了一种新的布置形式——变顶高尾水洞方案,当输水系统发生水力瞬变时,变顶高尾水洞尾水出流为明满交替流;采用无结构网格进行数值模拟分析,得到的结果将会更符合实际;另外,随着城市建设不断发展,供排水系统日趋复杂,当突降暴雨时,地面积水由于排水系统中存在气体而无法顺畅排出,水流流态由于地下系统气、液交界面的存在变得相当复杂,甚至会对地下排水系统的安全构成威胁.由此可见,很多涉及水、气两相交界,互相推动的水力瞬变问题亟待解决,而VOF方法作为处理自由面边界条件的重要方法,将在未来发挥越来越大的作用.
4 结论
VOF方法的发展和改进使两相、多相交界面的计算仿真越来越,流体输运方程由二维拓展至三维,差分格式复杂多样,边界适用条件从规则边界到不规则边界,计算网格从二维网格发展到三维网格、正交网格、非正交自适应网格、无结构网格等等.VOF法除提高自由面的模拟精度外,其处理的边界条件越来越复杂,适用范围越来越广,所解决的问题涉及化学、热能、机械、水利等众多学科和不同领域,在未来应用中将有更广阔的前景
VOF方法的发展和改进使两相、多相交界面的计算仿真越来越,流体输运方程由二维拓展至三维,差分格式复杂多样,边界适用条件从规则边界到不规则边界,计算网格从二维网格发展到三维网格、正交网格、非正交自适应网格、无结构网格等等.VOF法除提高自由面的模拟精度外,其处理的边界条件越来越复杂,适用范围越来越广,所解决的问题涉及化学、热能、机械、水利等众多学科和不同领域,在未来应用中将有更广阔的前景
